<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Figures 6.8-6.10: Quadratic smoothing</title>
<link rel="canonical" href="/Users/mcgrant/Repos/CVX/examples/cvxbook/Ch06_approx_fitting/html/smoothrec_cvx.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Figures 6.8-6.10: Quadratic smoothing</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 6.3.3</span>
<span class="comment">% Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX Argyris Zymnis - 10/2005</span>
<span class="comment">%</span>
<span class="comment">% Suppose we have a signal x, which does not vary too rapidly</span>
<span class="comment">% and that x is corrupted by some small, rapidly varying noise v,</span>
<span class="comment">% ie. x_cor = x + v. Then if we want to reconstruct x from x_cor</span>
<span class="comment">% we should solve (with x_hat as the parameter)</span>
<span class="comment">%        minimize ||x_hat - x_cor||_2 + lambda*phi_quad(x_hat)</span>
<span class="comment">%</span>
<span class="comment">% where phi_quad(x) = sum(x_(i+1)-x_i)^2 , for i = 1 to n-1.</span>
<span class="comment">% The parameter lambda controls the ''smoothness'' of x_hat.</span>
<span class="comment">%</span>
<span class="comment">% The first figure which is generated shows the original and</span>
<span class="comment">% the corrupted signals. The second figure shows the tradeoff curve</span>
<span class="comment">% obtained when varying lambda and the third figure shows three</span>
<span class="comment">% reconstructed signals.</span>
<span class="comment">%</span>
<span class="comment">% NOTE: This is not a good problem to use CVX on. By exploiting</span>
<span class="comment">% the sparsity in this case, we can solve this problem much more</span>
<span class="comment">% efficiently using least squares.</span>


randn(<span class="string">'state'</span>,0);

n = 4000;  t = (0:n-1)';
exact = 0.5*sin((2*pi/n)*t).*sin(0.01*t);
corrupt = exact + 0.05*randn(size(exact));

figure(1)
subplot(211)
plot(t,exact,<span class="string">'-'</span>);
axis([0 n -0.6 0.6])
title(<span class="string">'original signal'</span>);
ylabel(<span class="string">'ya'</span>);

subplot(212)
plot(t,corrupt,<span class="string">'-'</span>);
axis([0 n -0.6 0.6])
xlabel(<span class="string">'x'</span>);
ylabel(<span class="string">'yb'</span>);
title(<span class="string">'corrupted signal'</span>);
<span class="comment">%print -deps smoothrec_signals.eps % figure 6.8, page 313</span>

A = sparse(n-1,n);
A(:,1:n-1) = -speye(n-1,n-1);  A(:,2:n) = A(:,2:n)+speye(n-1,n-1);

<span class="comment">% tradeoff curve, figure 6.9, page 313</span>
nopts = 100;
lambdas = logspace(-10,10,nopts);

obj1 = [];  obj2 = [];

fprintf(<span class="string">'computing 100 points on tradeoff curve ... \n'</span>);

<span class="keyword">for</span> i=1:nopts

  lambda = lambdas(i);
  cvx_begin <span class="string">quiet</span>
    variable <span class="string">x(n)</span>
    minimize(norm(x-corrupt)+lambda*norm(x(2:n)-x(1:n-1)))
  cvx_end
  obj1 = [obj1, norm(full(A*x))];
  obj2 = [obj2, norm(full(x-corrupt))];

  fprintf(<span class="string">'tradeoff point %d\n'</span>,i);
<span class="keyword">end</span>;

figure(2)
plot(obj2,obj1,<span class="string">'-'</span>);  hold <span class="string">on</span>;
plot(0,norm(A*corrupt),<span class="string">'o'</span>);
plot(norm(corrupt),0,<span class="string">'o'</span>);  hold <span class="string">off</span>;
xlabel(<span class="string">'x'</span>);
ylabel(<span class="string">'y'</span>);
title(<span class="string">'||xhat-xcorr||_2 vs. ||D xhat||_2'</span>);
<span class="comment">%print -deps smoothrec_tradeoff.eps % figure 6.9, page 313</span>

<span class="comment">%three smooth signals, figure 6.10, page 314</span>
nopts = 3;
alphas = [8 3 1];
xrecon = [];

<span class="keyword">for</span> i=1:3
   fprintf(1,<span class="string">'Reconstructed Signals: %d of 3 \n'</span>,i)
   alpha = alphas(i);
   cvx_begin <span class="string">quiet</span>
    variable <span class="string">x(n)</span>
    minimize(norm(x(2:n)-x(1:n-1)))
    subject <span class="string">to</span>
        norm(x-corrupt) &lt;= alpha;
   cvx_end
   xrecon = [xrecon, x];

<span class="keyword">end</span>

figure(3)
subplot(311), plot(xrecon(:,1));
axis([0 n -0.6 0.6])
ylabel(<span class="string">'ya'</span>);
title(<span class="string">'||xhat-xcorr||_2=8'</span>);
subplot(312), plot(xrecon(:,2));
axis([0 n -0.6 0.6])
ylabel(<span class="string">'yb'</span>);
title(<span class="string">'||xhat-xcorr||_2=3'</span>);
subplot(313), plot(xrecon(:,3));
axis([0 n -0.6 0.6])
xlabel(<span class="string">'x'</span>);
ylabel(<span class="string">'yc'</span>);
title(<span class="string">'||xhat-xcorr||_2=1'</span>);
<span class="comment">%print -deps smoothrec_results.eps % figure 6.10, page 314</span>
</pre>
<a id="output"></a>
<pre class="codeoutput">
computing 100 points on tradeoff curve ... 
tradeoff point 1
tradeoff point 2
tradeoff point 3
tradeoff point 4
tradeoff point 5
tradeoff point 6
tradeoff point 7
tradeoff point 8
tradeoff point 9
tradeoff point 10
tradeoff point 11
tradeoff point 12
tradeoff point 13
tradeoff point 14
tradeoff point 15
tradeoff point 16
tradeoff point 17
tradeoff point 18
tradeoff point 19
tradeoff point 20
tradeoff point 21
tradeoff point 22
tradeoff point 23
tradeoff point 24
tradeoff point 25
tradeoff point 26
tradeoff point 27
tradeoff point 28
tradeoff point 29
tradeoff point 30
tradeoff point 31
tradeoff point 32
tradeoff point 33
tradeoff point 34
tradeoff point 35
tradeoff point 36
tradeoff point 37
tradeoff point 38
tradeoff point 39
tradeoff point 40
tradeoff point 41
tradeoff point 42
tradeoff point 43
tradeoff point 44
tradeoff point 45
tradeoff point 46
tradeoff point 47
tradeoff point 48
tradeoff point 49
tradeoff point 50
tradeoff point 51
tradeoff point 52
tradeoff point 53
tradeoff point 54
tradeoff point 55
tradeoff point 56
tradeoff point 57
tradeoff point 58
tradeoff point 59
tradeoff point 60
tradeoff point 61
tradeoff point 62
tradeoff point 63
tradeoff point 64
tradeoff point 65
tradeoff point 66
tradeoff point 67
tradeoff point 68
tradeoff point 69
tradeoff point 70
tradeoff point 71
tradeoff point 72
tradeoff point 73
tradeoff point 74
tradeoff point 75
tradeoff point 76
tradeoff point 77
tradeoff point 78
tradeoff point 79
tradeoff point 80
tradeoff point 81
tradeoff point 82
tradeoff point 83
tradeoff point 84
tradeoff point 85
tradeoff point 86
tradeoff point 87
tradeoff point 88
tradeoff point 89
tradeoff point 90
tradeoff point 91
tradeoff point 92
tradeoff point 93
tradeoff point 94
tradeoff point 95
tradeoff point 96
tradeoff point 97
tradeoff point 98
tradeoff point 99
tradeoff point 100
Reconstructed Signals: 1 of 3 
Reconstructed Signals: 2 of 3 
Reconstructed Signals: 3 of 3 
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="smoothrec_cvx__01.png" alt=""> <img src="smoothrec_cvx__02.png" alt=""> <img src="smoothrec_cvx__03.png" alt=""> 
</div>
</div>
</body>
</html>